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Abstract 

An approach to homogenization of high porosity metallic foams is explored. The emphasis is on the Alporas® foam 
and its representation by means of two-dimensional wire-frame models. The guaranteed upper and lower bounds on 
the effective stiffness coefficients are derived by the first-order homogenization with the uniform and minimal kinematic 
boundary conditions at heart. This is combined with the method of Wang tilings to generate sufficiently large material 
samples along with their finite element discretization. The obtained results are compared to experimental and numerical 
data available in literature and the suitability of the two-dimensional setting itself is discussed. 
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1. Introduction 

Highly porous metallic foams possess an extensive ap¬ 
plication potential. These materials feature high energy 
absorption, strength, and stiffness at very low weight, 
which makes them appealing for the automotive, aircraft, 
and defense industry, to name a few eh m Hi- In order 
to foster new application areas, a qualified understand¬ 
ing of the foam behavior and relevant predictive tools are 
required. Computer models, in particular, are regarded 
as a key ingredient to optimization of either/both the 
microstructure geometry or/and final products made of 
these materials [4]. A consensus regarding the approxi¬ 
mation of the behavior of porous solids seems to exist. 
Open-cell foam^] are usually represented with the three- 
dimensional beam models, while their closed-cell counter¬ 
parts require an addition of membrane elements acting as 
the cell walls [5]. Nevertheless, in the case of very thin 
walls even the behavior of closed-cell foams can be approx¬ 
imated with beams mm®- Following this assumption, 
Ashby and Gibson presented a three-dimensional beam 
model of the unit cell and derived solutions to the overall 
thermo-mechanical parameters based purely on material 
porosity [5]. Our goal is to explore the adequacy of a two- 
dimensional wired model discretized with beam elements. 
We build on the recent outcomes by Nemecek et al. [8], 
who developed a planar beam model based on the approx¬ 
imation of foam ligament geometry by the Voronoi tessel¬ 
lation. The objectives of the paper are threefold: 


^The author’s post-print manuscript of the article published in 
Computers & Structures , DOI: 10.1016/j.compstruc.2016.01.003 
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■As Banhart ,2j pointed out, sponge would be an appropriate 
term for open-cell highly porous solids. 


i) to investigate the influence of real geometry of foam 
ligaments and the proper size of computational mod¬ 
els [8] ; 

ii) to provide the upper and lower bounds on the effec¬ 
tive stiffness coefficients for the two-dimensional beam 
model via the first-order homogenization procedure; 

iii) to validate obtained results against the experimental 
and numerical results from [8’ and to question the as¬ 
sumption of the vanishing Poisson ratio made therein. 

The first objective is addressed with the help of Wang 
tilings, a concept recently introduced to Materials Engi¬ 
neering community |9]. In particular, the foam microstruc¬ 
ture is compressed within a set of smaller domains, called 
Wang tiles. The morphology of the tiles is designed such 
that the tiles are microstructurally compatible across the 
corresponding edges. As a result the reconstructed mate¬ 
rial microstructure remains continuous across the gridline 
of the regular lattice to which the tiles are accommodated 
during synthesis of a computational model, see Fig. J3^c). 
The method is extremely efficient in producing arbitrar¬ 
ily diverse ensembles of arbitrarily sized and geometrically 
consistent microstructure realizations in a fully stochastic 
setting M- Moreover, creating the finite element mesh 
on the level of tiles avoids expensive mesh generation of 
each microstructure realization. Altogether, we are able to 
reach for a proper size of the computational model, which 
is expected to be relatively large because of the infinite 
contrast in constituent properties m • 

As regards the second objective, we reformulate the first- 
order homogenization procedure for the wire-frame finite 
element models by means of macroscopic degrees of free¬ 
dom m The upper bound on the apparent stiffness is 
obtained from the ensemble of microstructures exposed to 
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the kinematic uniform boundary conditions (KUBC). The 
lower bound is rendered by applying the minimal kine¬ 
matic boundary conditions (mKBC) [13]. In order to com¬ 
pare our results with those reported in HI [2] [14j [8] , the 
bulk and shear moduli are derived from the homogenized 
stiffness matrices by assuming material isotropy. 

The paper is structured as follows. In the next chap¬ 
ter Alporas® aluminum foam, the material of interest, is 
characterized in brief. Section [3] covers the adopted mod¬ 
eling strategy addressing the first objective. In particular, 
Wang tilings are presented first together with their use 
in the context of numerical homogenization. Description 
of microstructure discretization into planar beams is also 
provided in this section. The first order homogenization is 
reviewed in Section [4] with emphasis on the beam model 
and the boundary conditions connected to the upper and 
lower bounds on the effective stiffness. Numerical results 
are provided in Section [5] and discussed in the last section. 

We use the following nomenclature within the text. 
Scalar quantities are denoted by plain letters, e.g. a or 
A , matrices by bold sans-serif font, e.g. a or A, and ten- 
sorial quantities by bold serif letters, e.g. a or A. In ad¬ 
dition, we adopt the standard Voigt matrix representation 
of symmetric second- and fourth-order tensors. 


In [8] the Mori-Tanaka homogenization scheme was em¬ 
ployed to estimate the values of effective Young’s modulus 
of 70 GPa and Poisson’s ratio of 0.35 of the bulk material. 
In order to obtain comparable data with [8], the values 
were taken here as the reference parameters of the solid 
(ligament) phase for the homogenization at the level of 
foam. 



Figure 1: Alporas® foam: a) Visualization of raw fi-CT data, b) 
image from optical microscope, c) binary image of polished cross- 
section, courtesy of Ondrej Jirousek 16] (a) and Jin Nemecek 0 
(b,c). 


3. Modeling strategy 


2. Alporas 

Alporas® is the closed-cell aluminum foam commer¬ 
cially produced by Shinko Wire Company, Ltd. m with 
porosity between 90 % and 92 % pQ. The pore space 
is composed of approximately polyhedral voids with the 
mean size of around 3 mm. The main field of its appli¬ 
cation lies in the energy absorption structures. Predomi¬ 
nantly, Alporas® is used in production of self-supporting 
sound absorbers m- The sandwich plates comprising of 
Alporas® foam core were also reported to outperform com¬ 
mon steel plates under blast loading in e.g. ng. 

Regarding the linear elastic behavior, the design guide to 
metallic foams [I] suggests the overall parameters enlisted 
in Tab. [l] However, the region of linear elastic behavior 

Table 1: Properties of Alporas® foam 1 . Here K, G and E stand 
for the bulk, shear and Young moduli, respectively, and v is the 
Poisson ratio._ 


K [GPa] 

G [GPa] 

E [GPa] v [-] 

0.9 - 1.2 

0.30 - 0.35 

0.4 - 1.0 0.31 - 0.34 


is rather limited. Nonlinear or inelastic behavior develops 
soon after the loading onset due to wall buckling, contacts 
between ligaments, and local yielding, which makes inter¬ 
pretation of overall elastic parameters from experimental 
data a challenging task. 

Two distinct phases, (i) the aluminum rich phase, and 
(ii) the phase of precipitates consisting of titanium and 
calcium, were identified by Nemecek et al. at the level of 
cell walls for the investigated specimen displayed in Fig.[^ 


Numerical homogenization is closely connected with the 
notion of Representative Volume Element (RVE), a char¬ 
acteristic sample of a microstructure large enough to pro¬ 
vide us with the effective material properties. From this 
viewpoint, a method capable of rendering stochastic mi¬ 
crostructure realizations of arbitrary sizes is desirable. 

Our approach rests on Wang tiling, a multicell general¬ 
ization of the Periodic Unit Cell concept [9], see Fig. [2j We 
adopt the premise that a microstructure represented to a 
high degree of accuracy automatically provides us with a 
broad range of effective properties mm- Starting with a 
binary two-dimensional image of a material, we compress 
the microstructure into a set of Wang tiles. Subsequently, 
the Finite Element (FE) mesh composed of planar beams 
conforming through the edges of identical codes is gener¬ 
ated within each tile. Finally, in the numerical homoge¬ 
nization we benefit from the ability of the tiling concept to 
produce ensembles of microstructure realizations and cor¬ 
responding computational models of arbitrary sizes using 
a stochastic tiling algorithm. 

3.1. Wang tilings 

Wang tiles are square domino-like pieces with codes as¬ 
signed to their edges, Fig.[3ja). The codes represent bind¬ 
ing constrains such that only the tiles with the same codes 
on the adjacent edges may be placed next to each other 
when assembled into a tiling. By definition, the tiles can 
be neither rotated nor reflected m , hence the two tiles of 
identical code sequences mutually rotated by f are consid¬ 
ered as different. The term tile set denotes the union of all 
tiles available for the tiling synthesis. The set is character¬ 
ized with the number of tiles and the number of distinct 
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Figure 2: Illustration of adopted modeling strategy. 


codes on their horizontal and vertical edges. The compo¬ 
sition of codes within the tile set governs the probability 
of occurrence of each tile in the tiling. 

Specific families of tile sets can be identified. Each of 
the families usually comes with a unique tiling algorithm. 
A great attention, especially from Discrete Mathematics 
community, has been paid to the family of aperiodic sets 
in the last 50 years. However, in agreement with Com¬ 
puter Graphics applications [20], we prefer stochastic sets 
to strictly aperiodic ones as the former provide us with 
more freedom in their design. The stochastic sets also cor¬ 
respond better with the purpose of modeling of materials 
with random/disordered microstructure. 

Tiles of a stochastic tile set can be easily assembled into 
tilings by making use of the stochastic tiling algorithm 
proposed by Cohen et al. [20] . First, a regular lattice of 
desired size is created, see Fig. The grid is sequen¬ 
tially filled up with tiles in either column-by-column or 
row-by-row order. At each step, a tile is randomly chosen 
from a subset of tiles that satisfy edge constrains given by 
the previously placed tiles and positioned. The random¬ 
ness is guaranteed by the existence of at least two distinct 
tiles for each admissible edge code combination on the top 
and left-hand side edges. A schematic of the procedure 
is displayed in Fig. [3](c) . The tiles 4 and 3 from the tile 
set, Fig. §b), form the unique subset for one step in the 
stochastic tiling algorithm sketched in Fig. [3][c). Either of 
the tiles is randomly chosen and fixed at the spot indicated 
with the asterisk mark. 
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Figure 3: (Color online) Illustration of: a) Wang tile, b) set of 8 
Wang tiles, c) single step of stochastic tiling algorithm. 


3.2. Tile-based micro structure compression 


In order to make the Wang tile concept applicable to 
modeling of random heterogeneous materials, microstruc- 
tural information has to be attributed to each tile such 
that a tiled domain (a tiling) corresponds to the original 
microstructure and the microstructure continuity is not 
violated across the congruent edges. By analogy to the 
Periodic Unit Cell (PUC) morphology design [18J [21] . a 
way to meet those requirements is to employ optimiza¬ 
tion methods, e.g. Simulated Annealing [9]. Although 
this approach is fully general, the computational cost is 
tremendous especially for larger sets, high microstructure 
resolution, multiphase media, and/or multi-point spatial 
descriptors. Therefore, in this study we employ the auto¬ 
matic design procedure proposed originally for compres¬ 
sion of digital textures [20] and explored later in materials 
modeling m ■ 

In brief, we take a number of samples from the refer¬ 
ence specimen of the investigated microstructure. Each of 
the samples is associated with a unique edge code from 
the tile set. Each tile is then constructed as a square 
cut, rotated by f, out of partially overlapping samples 
that are positioned according to the tile edge codes they 
are associated with BOj. Since the cut goes diagonally 
through the samples, the microstructure continuity is au¬ 
tomatically preserved as the sample leftovers are used as 
microstructural information for the tiles having the same 
codes on the opposite edges. The samples are seamlessly 
fused in the tile interiors by making use of the Image Quilt¬ 
ing algorithm due to Efros et al. [22]. We have augmented 
the procedure by optimizing the design input parameters 
to achieve maximum proximity in spatial statistics of syn¬ 
thesized tilings and the original microstructure. Namely, 
volume fraction, two-point probability and two-point clus¬ 
ter functions were considered. This approach appears ex¬ 
tremely efficient from the viewpoint of both the quality of 
compressed microstructures and computational overhead, 
see [10] for details. Thus, after the compression based on 
the automatic tile design algorithm we have a binary rep¬ 
resentation of the Alporas microstructure in the form of a 
set of Wang tiles at hand. 
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3.3. Discretization 

Validity of three-dimensional beam models for predict¬ 
ing elastic behavior of foams has been reported several 
times, e.g. B ID 13- In the recent work [8], Nemecek 
et al. presented promising results also for the two- 
dimensional beam representation of foam motivating our 
modeling strategy. In their approach the foam microstruc¬ 
ture was represented as the Voronoi wired model consist¬ 
ing of straight beams with the centroids of individual pores 
taken as the Voronoi diagram seeds. The wired model was 
preferred to the two-dimensional finite element analysis as 
the latter would have required modifications to ligament 
walls thickness in order to fit the volume fraction from 
the weighting experiment. Its difference against the vol¬ 
ume fraction determined from the area of ligaments in the 
reference planar scan was attributed to the preparation 
procedure, in which the saw cut was not perpendicular to 
most of the cell walls yielding the apparent solid phase 
volume be higher than the true one. 

The size of the model used by authors of [8] was limited 
by the extent of the given scanned image. However, this 
is not the case of Wang tiles, which enable to investigate 
series of models yielding adequate RVE dimensions. In ad¬ 
dition, the beam geometry conforming with the real shape 
of ligaments allows us to explore the inaccuracy result¬ 
ing from the Voronoi tessellation based idealization. We 
started with the foam geometry compressed within the set 
of 32 tiles with four distinct edge codes on vertical and hor¬ 
izontal edges (8 in total), Fig. 0- First, a parametric beam 
discretization of tile geometries was created manually by 
employing a specifically designed GUI tool written in Mat- 
lab environment^ Fig. S>- It enabled us to control the 
mesh compatibility across tile edges and to fix ligaments 
disrupted by quilting or faults due to the sample prepara¬ 
tion, Fig. [I] Computational domains were assembled from 
the individual tiles by making use of the stochastic tiling 
algorithm described in Section |3.1[ the nodes on the coin¬ 
cident edges were identified, merged, and the entire record 
renumbered to maintain the sparsity of algebraic systems. 

The cross-sectional parameters of individual beams were 
determined by analogy to [8]. Beams were assumed 
straight and prismatic with a rectangular cross-section of 
the unit width. The total length of beams L e is given by 
the mesh geometry. Hence, the only degree of freedom, the 
depth of the cross-section iL, is governed by the volume 
fraction of the solid phase p and RVE volume/area \Q\. It 
reads as 


H = 




(i) 


Together with the assumption of unit beam width, Eq. 0 
further yields expressions for the cross-sectional area A 
and the second moment of area / 


A = H, 



( 2 ) 


2 available at http://mech.fsv.cvut.cz/~doskar/software/ 
MeshMaker.zip 


In order to analyze precisely the effect of beam geome¬ 
try, a separate beam model, Fig. [5^, was created directly 
from the reference geometry (referred to as real geometry 
mesh onwards) in our in-house developed GUI tool. It was 
accompanied also with the microstructure represented as 
the Voronoi diagram, Fig. 0, to compare with the results 
published in [8]. For visual comparison, a tiling-based ge¬ 
ometry of similar size (composed of 6 x 6 tiles) is shown 
in Fig. [5]i. 


4. Numerical homogenization 


In general, the homogenization process provides a homo¬ 
geneous substitute for a heterogeneous material composed 
of constituents of different properties. Limiting the expo¬ 
sition to linear elasticity, with Hooke’s Law 

cr(x) = D(as) : e{x) , for x G U , (3) 

relating local stresses cr{x) to strains s(x) and U being the 
microscale domain of interest representing a macroscopic 
material point, we wish to replace the heterogeneous ma¬ 
terial stiffness tensor D(cc) with its homogeneous counter¬ 
part D hom independent of microscale coordinates x. An 
equation analogical to Eq. 0, namely, 

£ = D hom : E , (4) 


then relates the macroscopic strains E and stresses X, de¬ 
fined as spatial averages (•) 


E = (e(x)} = T“~7 f e(x) dx and 

bd Jn 

£ = -Tfiif <T ( a: ) da3 ’ 


(5) 


respectively. 

Homogenization of parameters describing linear prob¬ 
lems is a well resolved area nowadays. For simple mi- 
crostructural compositions and relatively low contrast in 
constituent properties analytical homogenization schemes 
yield accurate estimates, often based on the Eshelby solu¬ 
tion, or narrow bounds, e.q. Voigt-Reuss-Hill or Hashin- 
Shtrikman, of the sought homogenized property [23], and 
references therein]. On the contrary, for materials with 
high contrast in properties of its constituents the actual 
microstructural geometry, simplified in the analytical ho¬ 
mogenization methods, influences significantly the mate¬ 
rial response. Consequently, the analytical methods are 
suitable only for quick precursory estimates of the over¬ 
all material stiffness. In the particular example of highly 
porous foams, even the closest estimates obtained with the 
Mori-Tanaka and the differential schemes differ from the 
experimental data by 20 % [24j 23]. 

Hence, to determine the overall properties we resort to 
a numerical homogenization that, albeit more computa¬ 
tionally intensive, gives us an insight into the size of the 
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Figure 4: Alporas® microstructure compressed within 32 tiles: (a) microstructure of individual tiles, (b) corresponding wired model (beam 
discretization of ligaments). 



Figure 5: (Color online) Comparison of various representations of Alporas® geometry, a) reference microstructure (courtesy of Jiff 
Nemecek [8]), b) Voronoi mesh, c) mesh of real geometry, d) tiling-based representation of beam mesh (6x6 tiling). The red line high¬ 
lights boundary beams such that they coincide with dQ. 


related RVE and mechanics of deformation of the inves¬ 
tigated microstructure. In particular, we adopt the first- 
order numerical homogenization method described com¬ 
prehensively in work of Michel et al. P3- First, the gen¬ 
eral concept is recalled in the tensorial notation and the 
variational setting. Then, we reformulate the approach for 
a two-dimensional beam model, using the matrix notation. 
Numerical homogenization, unlike the analytical schemes, 
calls for a suitable computational domain - RVE, and for 
specific boundary conditions in order to render bounds of 
the effective property. We address these issues in Sec¬ 
tion |4.3| and discuss incorporation of different boundary 
conditions by means of the Lagrange multipliers. 

4-1- First order homogenization 

In the first order homogenization, the total displacement 
field u(x) can be decomposed as m 

u(x) = u E {x) + u*(x) = E • x + u*(x ), force ell, (6) 

where u E (x) is coupled with the homogeneous strains E 
and u*(x) are displacement fluctuations induced by pres¬ 
ence of material heterogeneity. Applying the symmetric 
gradient operator V s = \ (V + V T ) to the displacement 
fields in Eq. ([6| provides us with the total strains in the 


form 

e(x) = V s u(x ) = E + V s u*(x ), forxeQ. (7) 

Thus, the strain energy density on microscale reads 

W{x,E, V s u*(x)) = 

1 ( 8 ) 

- (E + V s u*{x)) : DO): (E + V*u*(x )). 

For the purpose of variational formulation, we distin¬ 
guish between the true fields E and u* and the corre¬ 
sponding test fields denoted with tilde. Let C be a set 
of all admissible fluctuation displacement fields complying 
with the boundary conditions, further specified in |4.3| For 
arbitrary test fields u* G C and E E where d is the 

dimension of the model under consideration (d = 2 in our 
particular case), we define the macroscopic strain energy 
density W as 

W{E, u*) = E, V s u*(x)4 ■ (9) 

Assume the stress 5] to be prescribed on macroscale. The 
elastic energy potential J(E,u*) for a macroscopic mate- 





















rial point is then specified in the form 


macroscopic strain energy can be written as 


J(E,u*) = W(E,u*)-E: E. 


( 10 ) 


The true fields E and u * are obtained by minimizing 
the latter equation following the minimum total potential 
energy principle. Finally, combining the stationary condi¬ 
tions of Eq. (10) and Eq. 0 results in the expression 


, d 2 W 

D hom = t E u *\ 

dE 2 


(ii) 


4-2. Discretization - Macroscopic Degrees of Freedom 

As discussed in the previous sections, we assume D to be 
a planar model discretized by straight beams. The model 
is embedded in the x-z plane, thus the unknowns of the 
beam model are nodal displacements u , w along x, z axes, 
and rotations (j) about y. The dual quantities are normal 
forces TV, shear forces V and bending moments M acting 
at beam ends (denoted by subscript indices 1 and 2). It 
holds 

fe = K e u e , (12) 

where f e = {TVi, Vi, M 1 , TV 2 , M 2 } T , K e is the element 
stiffness matrix calculated according to the Timoshenko 
theory [25] . and u e = {ui, Wi, 0i, u 2 , w 2 , (/) 2 } T . 

The displacement decomposition in Eq. 0, rewritten in 
the matrix form, gives 


Ue = [I 



(13) 


where I is the identity matrix of 6 x 6 entries and 


A 


T 

e 


x\ 0 0 X 2 0 0 

0 zi 0 0 Z 2 0 

\z\ \x\ 0 \z 2 \x 2 0 


(14) 


couples the nodal degrees of freedom (DOFs) with the ap¬ 
plied macroscopic strains E = {E x , E Zl r xz } J . Thus, the 
components of the macroscopic strain play the role of ad¬ 
ditional macroscopic degrees of freedom of the discretized 
system, cf. [la¬ 
in the variational setting, contributions of individual fi¬ 
nite elements to the overall strain energy density are 


W e ( E,u*) = 



Kjx 


(15) 


where K® x 
form 


is the extended element stiffness matrix of the 


K 


ex 

e 


K e K e A e 
AjK e AjK e A e 


(16) 


With the extended discrete 


system (15) at hand, the 




m 1 


where u* = Agf. x u*, K ex takes the form 


K ex = A Kf = 

e=l 


K ex 

■Ml 

K ex 

■^21 


K ex 

^12 

K ex 

^22 


( 17 ) 


(18) 


and is the finite element assembly operator. The 

submatrices K® x arise from the assembly of the submatri¬ 


ces on the right-hand side of Eq. (16) and n e denotes the 
number of elements. 

4-3. RVE size and boundary conditions 

The above procedure yields the homogenized parame¬ 
ters providing the computational domain coincides with 
the RVE related to the investigated microstructure and 
physical phenomenon. According to Hill [26] . RVE is the 
sample of a material that incorporates enough geometri¬ 
cal information to render the effective stiffness parameters 
regardless the prescribed boundary conditions as long as 
they would result in the uniform stress and strain states 
if the microstructure were homogeneous. For smaller do¬ 
mains, here called Statistical Volume Elements (SVEs), 
Eq. (11) gives only the boundary condition dependent ap¬ 


parent stiffnesses. Thus, when homogenizing, the size of U 
is usually increased step by step until the apparent prop¬ 
erties resulting from different boundary conditions either 
coincide, i.e. meeting Hill’s criterion [2S], or the variance 
in the sought quantity computed from different realizations 
of the same size reaches a given limit [27] . 

Among the admissible boundary conditions yielding the 
uniform strain in the homogeneous medium, three are com¬ 
monly used. The periodic boundary conditions show the 
fastest convergence as the size of SVEs approaches that 
of RVE, however they render only the estimate of the ef¬ 
fective stiffness and are cumbersome to impose if the mi¬ 
crostructure is in general non-periodic, which is also the 
case of Alporas®. The remaining two, Kinematic Uniform 
Boundary Conditions (KUBC) and Static Uniform Bound¬ 
ary Conditions (SUBC), are known to result in the upper 
and lower bounds of the apparent stiffness coefficients of 
each SVE, refining so the Voigt-Reuss-Hill bounds [271128]. 
To the best of our knowledge, the exact form of SUBC is 
unknown for beam models. We thus circumvent this draw¬ 
back prescribing the Minimal Kinematic Boundary Con¬ 
ditions (mKBC) due to Mesarovic and Padbidri [13!, since 
these were shown to equal SUBC [29] . 

To preserve the structure of the resulting algebraic sys¬ 
tem, both boundary conditions are prescribed as an addi¬ 
tional constraint^ in the form 


CQ* = 0 


(19) 


3 We assume that the assembly operator A^ 1 does not involve 
boundary conditions. 
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defining so the kinematically admissible nodal displace¬ 
ments and rotations stored in u*. The boundary dfl has 
been determined from the beams fully contained in the 
domain, highlighted in red (light gray in BW version) in 
Fig. [5] 

In the case of KUBC, the displacement fluctuations have 
to vanish at the boundary, i.e. u*(x) = 0 for x G dfl . 
In the discrete beam model we prescribe zero fluctuation 
displacements u and v for all nodes at the boundary, the 
corresponding rotations <f> are left unknown. With n dQ 
boundary nodes and n tot nodes in total, the matrix of con¬ 
straints C KU for the KUBC case consists of [2 n dn x 3n tot ] 
components and is defined as 

C^ = I 1’ forties, (20) 

[ 0, otherwise, v ' 

where Q is the set of all code numbers corresponding to 
the fixed fluctuation displacements. 


In contrast to KUBC, mKBC [29j do not require the 
fluctuation field to vanish at the boundary point-wise but 
they enforce the field to vanish on average. The definition 
of the macroscopic strain tensor E , Eq. JT along with the 
local strain field decomposition, Eq. 0. yields the mini¬ 
mal condition posed on u* (also called the normalization 
condition [30]) in the form 

0= H/V'“' d * = H/ M "®' fi ' dx ' (21) 


where n is the outer normal and (g) s denotes the sym¬ 
metric part of the tensor product defined as (n 0 s u*) = 
\ (n (8) u* + u* (8) n) . In addition to the three conditions 
defined by Eq. ( |2Tj ) (for the two-dimensional setting), we 
prescribe zero fluctuations u* and w* at the bottom-left 
and w* at the bottom-right corners to restrain the rigid 
body modes. Plugging the shape functions of the beam 
element in Eq. © along with the three additional condi¬ 
tions mentioned above give rise to the matrix of cons traints 
C mK with [(3 + 3) x 3 n dn ] entries, see|~^ 
ditional details. 


Appendix A 


for ad- 


form the following algebraic system for the true fields 



( 24 ) 


Finally, condensing u* and A from Eq. (24) yields the ho¬ 
mogenized stiffness matrix D hom equivalently to Eq. (11). 
For large and sparse systems of linear equations direct con¬ 
densation of Eq. (24) would be inefficient. We thus com¬ 
pute individual columns of D hom in three consecutive steps 
as the responses ZW to the unit load cases defined as 
EW = {1,0,0} T , E( 2 ) = {0,1,0} T , and E< 3 ) = {0,0,l} T . 


4-4- Extraction of isotropic parameters 


The approach outlined above leads to the homogenized 
material stiffness matrix with components encoding the 
anisotropic behavior. If the degree of anisotropy is small, 
it is convenient to describe the constitutive behavior of 
the homogenized material with scalar quantities K and 
G, see Tab. [l] To this goal, we perform spectral analysis 
of D hom to determine the effective moduli from its eigen¬ 
values m- For the plane strain conditions, the stiffness 
matrix of isotropic material D lso and its eigenvalues A lso 
in terms of K and G read 


D iso = 


~K + 

k-\g 

K - |G 
K+lG 

0" 

0 

, A iso = J 

0 

0 

G 

1 


G 1 

2 G } 

|G + 2K J 


(25) 


Note, that the min —max ordering of the A lso components 
holds for non-auxetic materials only, i.e. for those with 
v > 0. Since the calculated overall stiffness matrix D hom 
and the matrix eigenvalues A hom may not obey Eq. (25) 
exactly we employ the Least Square Method to obtain the 
moduli as 


J G hom \ 

" 1 

5 

1 ol 

(A' 110111 ] - 

i 

4 1 

30 2. 


see |Appendix B| for additional details. 


(26) 


5. Results 


Having the matrices C KU and C mK at hand, the con¬ 
strained minimization of Eq. (10) can be reformulated by 
making use of the Lagrange multipliers A, which act as re¬ 
actions and generalized reactions for KUBC and mKBC, 
respectively. To this goal, we define the energy potential 
J(u*,E,A) similarly to Eq. (10) 


J(u*,E,A) 



The stationary conditions posed on J, namely, 


d 

<9u* 



d -\ 

4 F = 0 , 

<9E lE=E 


A J| -0, (23) 
8X lx=x 


The cross-sectional characteristics of beams, Tab. [2j for 
the three different geometries, Fig.[5j were derived through 
the approach outlined in Section [3] with the target liga¬ 
ment volume fraction set to 8.6 % [8]. Note that the depth 
of beams (corresponding to the cross-section area A , re¬ 
call Eq. 0 ) in the tiling-based representation is about 
10 % smaller than for the previous two. We attribute this 
reduction to the automatic tile design, namely local mod¬ 
ifications to the geometry within the overlap region, and 
corrections induced by manual mesh generation which may 
resulted in longer beams and the corresponding smaller 
depth H. 

We have computed the apparent parameters for square 
computational domains/tilings of sizes ranging from lxl 
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Table 2: Cross-sectional characteristics of beams conforming with 
weighted volume fraction. 



Voronoi 

real-geometry 

tiling-based 

A [m 2 ] 

1.431 x 10“ 04 

1.469 x 10“ 04 

1.317 x 1CT 04 

I [m 4 ] 

2.443 x 10“ 13 

2.643 x 10“ 13 

1.904 x 10“ 13 



Figure 7: (Color online) Coefficient of variation of the upper and 
lower estimates on K and G. 

up to 300 x 300 tiles. For each size, we generated 32 differ¬ 
ent stochastic realizations so that the number of realiza¬ 
tions equals the number of distinct tiles in the set. Thus, 
results for 1 x 1 tilings correspond to the averages over the 
individual tiles. 

The means of the effective stiffness coefficients for mod¬ 
els consisting of 300 x 300 tiles are enlisted in Tab. [^to¬ 
gether with the homogenized parameters arising from the 
Voronoi and real-geometry meshes as presented in Fig. [5] 
As expected, the upper bounds generated with KUBC and 


Table 3: Apparent elastic parameters of individual models. 




J^hom 

^hom 

^hom 

^hom 



[MPa] 

[MPa] 

[MPa] 

H 

Voronoi 

KUBC 

1435 

238 

81 

0.472 


mKBC 

1027 

51 

17 

0.492 

real-geometry 

KUBC 

203 

220 

84 

0.319 

mKBC 

61 

34 

12 

0.409 

tiling-based 

KUBC 

107(144) 

42(57) 

15(20) 

0.435(0.434) 

mKBC 

104(141) 

38(53) 

13(18) 

0.439(0.438) 


the lower bounds rendered by means of minKBC converge 
to each other with increasing size of computational do¬ 
mains, Fig. [6] even for tilings of 300 x 300 tiles they do 
not fully coincide. Nevertheless, we consider the difference 
between the bounds acceptable. This conclusion can be 
justified also in terms of the variance of the bounds. Fig. [7] 
shows the asymptotic decay of coefficient of variation with 
the values less than 1 %o in the vicinity of 300 x 300 tiling. 

Since we failed to predict the reported stiffness param¬ 


eters of Alporas® foam, Tab. [I] with the cross-sectional 
characteristics derived for tiling-based geometries, we per¬ 
formed the same computations also for cross-sectional pa¬ 
rameters derived for the real-geometry mesh, Tab. [2] The 
results, in brackets in Tab. [3j were approximately 40 % 
higher reflecting nearly linearly the 40% increase in the 
second moment of area. However, they do not fall within 
the reported range of overall Alporas® properties either. 

6. Conclusions and discussion 

We presented a modeling strategy for Alporas® foam 
having at heart the synthesis of stochastic microstructure 
realizations based on the concept of Wang tiling. Besides 
the standard upper bounds on the effective stiffness given 
by Kinematic Uniform Boundary Conditions, we obtained 
the guaranteed lower bounds by means of minimal kine¬ 
matic boundary conditions. With this formulation, shown 
to equal the common Static Uniform Boundary Conditions 
in [29], we can avoid the question of an appropriate beam 
loading that would yield macroscopically uniform stress 
field. Hence, both bounds were rendered prescribing rele¬ 
vant boundary displacement. 

Based on data in Fig. [6] we conclude that RVE for the 
two-dimensional model of highly porous materials should 
be about hundreds of the characteristic pore diameter 
length in dimensions. This finding corresponds with state¬ 
ments made in m regarding the minimal size of RVE 
in the case of infinite contrast of phase properties. On 
the other hand, it contradicts the recommendation of 
Ashby et al. [1 who propose RVE size of approximately 
seven times the mean pore diameter. However, this rec¬ 
ommendation is given for three-dimensional samples and 
does not need to be valid for planar analyses. 

Comparison of the homogenized stiffness coefficients 
with the reference values reported by Ashby et al. [T], Tabs. 
[l]and[3] leads us to unambiguous conclusions. Despite the 
several times reported aptness of the spatial wired model 
n s m, it can be conjectured that the adopted planar 
beam representation has a limited capability in predicting 
the complex behavior of Alporas® foam. Possibly, it lacks 
the stiffness contribution from the out-of-plane beams and 
membranes as well as the cell walls parallel to the inves¬ 
tigated plane, though, the membrane contribution of cell 
walls was reported negligible in the case of high-porosity 
foams [7]. 

The qualitative analysis of the impact of the geome¬ 
try representation clearly shows that the Voronoi mesh 
adopted in |8j leads to the overestimated value of bulk 
modulus. Assuming only the volumetric deformation, the 
axial stiffness of straight beams dominates the behavior of 
the Voronoi model resulting in nearly incompressible be¬ 
havior, whilst in the case of real-geometry and tiling-based 
meshes the axial and bending stiffness contribute equally. 
This explains the results reported in |8j where the authors 
considered only the volumetric excitation. On top of that, 
based on their experimental observations, they assumed 
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(a) (b) 



(c) (d) 


Figure 6: (Color online) Relationship between homogenized elastic parameters and dimensions of computational domains: (a) Bulk modulus K , 
(b) Young’s modulus E , (c) Shear modulus G, (d) Poisson’s ratio v. 


negligible Poisson’s effect and equated Young’s modulus 
to the oedometric component of the homogenized stiffness 
matrix. However, if the homogenization procedure was 
performed carefully, they would arrive at similar results 
as in Tab. 02 It is worthwhile to note that zero Poisson’s 
ratio is in contradiction also with the characteristics listed 
in Tab. |TJ see JY for further details. 

Despite the incapability of the two-dimensional beam 
model to capture behavior of complex materials such as 
highly porous metallic foams, the concept of Wang tilings 
proved to be a valuable technique to generate large com¬ 
putational domains of desired geometrical characteristics. 
Moreover, it allows to circumvent meshing of every single 
microstructure realization if the dicretization mesh is con¬ 
structed on each tile separately such that it satisfies the 
given compatibility constrains. There is a promising pos¬ 
sibility to extend the presented modeling strategy to three 
dimensions via Wang cubes, a spatial variant of Wang tiles, 
along with the spatial Voronoi tessellation to compress the 
true geometry of Alporas® foam, and to perform a similar 
analysis using shell elements instead of prismatic beams. 
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Appendix A. Minimal kinematic boundary condi¬ 
tions for beam model 


In the two-dimensional wireframe setting, Eq. (21) 


defining Minimal Kinematic Boundary Conditions can be 
rewritten into the following form 


° = E [ 

Jlw 


(i) 


jnl^w(x?) + n^ > u{x t ) j 


W) 

w(x e ) 
«„ 


dx* 


(A.l) 


where nf n stands for the number of the boundary beams, 
is the length of the i-th boundary beam, n\^ and 
are the components of the outward unit normal vector nW, 
and x £ is the local coordinate measured along the beam 


axis. Eq. (A.l) allows to compute the i-th boundary beam 
contribution to the constraint matrix C mK , Eq. (19) as 


^mK,(i) 


n 


W 

l 

0 




W 


(i) ( i ) 

io J n) 


J £ ~ g M J g ~ £ 


(A.2) 
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where the matrices T^ _s and J g ~ £ transform the corre¬ 
sponding quantities between the local and the global co¬ 
ordinate systems. Namely, 


T^ s 


cos a W — sin a W 


V W 


V W 


(A-3) 


and 


T s_ ^ = 


cos aW sin aW 0 0 0 0 

— sinaW cosaW 0 0 0 0 

0 0 1 0 0 0 

0 0 0 cos ck W sinaW o 

0 0 0 — sinaW cosaW o 

0 0 0 0 0 1 


(A.4) 


where is the directional angle of the i-th beam. The 
matrix M 


M = 


0 

2 u 

0 ^ 

u 2 


0 ^ 0 

■^2 T (i) 

U 2 


12 


0 

—IT) 2 

12 


(A.5) 


The least square method applied to Eq. (B.l) results in 


minimizing the distance between the calculated eigenval¬ 
ues A hom of the homogenized stiffness matrix and A lso in 
the Euclidean norm 

|| A iso _ A hom||2 = ( A iso _ A homy ^iso _ Ahom ^ _ ( B-2 ) 

Plugging Eq. into Eq. ( |B.2[ ) and minimizing it with 

respect to d yields the parameters K hom and G hom explic¬ 
itly approximated as a projection of A hom onto the range 
of A 


d h om _ 1 A T A hom = 


f 0 

__4_ 1 

30 2 


A 


hom 


(B.3) 


where (A t A) 1 A T is the left pseudoinverse of A. 

Note that applying the above procedure in plane-stress 
conditions leads to non-linear relation between the eigen¬ 
values and the elastic parameters 


A" 


{G 2 G 


18 KG 1 1 
4G+3K J 


(B.4) 


is obtained by plugging the the shape functions for u and 
w, e.g. [32], 


As a result, no explicit formula can be derived and an 
iterative Least square method must be used instead. 


K 1 (0-i -O 
K 2 (0-e, 

N .T (0 = 1 

N^ 1 (£) = 

w 1 + 2 K 

K 2 (0 = r ^:[(2^ + 3C 2 -2^ 3 ], 

n t (0 = r ^[< + ( 1 -«)C 2 -4 3 ], 


1 T 2 Kj 
L 


[(1 + 2k) — 2k£ — 3£ 2 + 2£ 3 ], 

[-(i + K)e+(2 + K)^ 2 -e 3 ], (A - 6 ) 


where £ is the normalized coordinate £ = ^y, k = 
and k is the Timoshenko shear coefficient, in the approxi¬ 
mation of the local beam displacements 
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N U2 
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N Wl 
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N* 1 
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]\[W2 
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w i 
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uu 


(A.7) 


and evaluating the integral in Eq. (A.l). 


Appendix B. Least square method 

The eigenvalues of the isotropic plane-strain stiffness 
matrix written as a linear combination of parameters K 
and G read as 


A lso = Ad = 


- 2 
.3 z 


(T 


(B.l) 
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